#include "libmri.h"

#include <algorithm>
using std::min;
using std::max;
/*

void kreduce3d (cplxf *KKt, cplxf *K, int ks1, int ks2, int ks3, int n, int m)
{
	int kk1 = 2*ks1-1, kk2 = 2*ks2-1, kk3 = 2*ks3-1;

  #pragma omp parallel
  {
	#pragma omp for schedule(static)
	for (int i = 0; i < kk1*kk2*kk3*n*n; ++i){ KKt[i] = 0.0f; }

	#pragma omp for schedule(static)
	for (int i = 0; i < n; ++i){
	  for (int j = 0; j < n; ++j)
		for (int k = 0; k < m; ++k)
		{
			cplxf *KKtij = &KKt [(i + n*j)*kk1*kk2*kk3],
			        *Kik =   &K [(i + n*k)*ks1*ks2*ks3],
			        *Kjk =   &K [(j + n*k)*ks1*ks2*ks3];

			for (int i3 = 0; i3 < kk3; ++i3)
			  for (int i2 = 0; i2 < kk2; ++i2)
	 			for (int i1 = 0; i1 < kk1; ++i1)
				{
					int j3_0 = max(0, i3-ks3+1), j3_n = min(i3+1, ks3);
					int j2_0 = max(0, i2-ks2+1), j2_n = min(i2+1, ks2);
					int j1_0 = max(0, i1-ks1+1), j1_n = min(i1+1, ks1);

					for (int j3 = j3_0; j3 < j3_n; ++j3)
					  for (int j2 = j2_0; j2 < j2_n; ++j2)
						for (int j1 = j1_0; j1 < j1_n; ++j1)
						{
							KKtij [i1 + kk1*(i2 + kk2*i3)] += 
									Kik [(i1-j1) + ks1*
									    ((i2-j2) + ks2*
									     (i3-j3))] *
									conj (Kjk [(ks1-1-j1) + ks1*
									          ((ks2-1-j2) + ks2*
									           (ks3-1-j3))]);
						}
				}
		}
	}
  }
}
*/

void kreduce3d (cplxf *KKt, cplxf *K, int *ksize, int n, int m)
{
	int ks1 = ksize[0], ks2 = ksize[1], ks3 = ksize[2];
	int kk1 = 2*ks1-1, kk2 = 2*ks2-1, kk3 = 2*ks3-1;

	#pragma omp parallel
	{
		#pragma omp for
		for (int i = 0; i < kk1*kk2*kk3*n*n; ++i){ KKt[i] = 0.0f; }

		#pragma omp for
		for (int i = 0; i < n; ++i)
		{
		  for (int j = 0; j < n; ++j)
			for (int k = 0; k < m; ++k)
			{
				cplxf *KKtij = &KKt [(i + n*j)*kk1*kk2*kk3],
				        *Kik =   &K [(i + n*k)*ks1*ks2*ks3],
				        *Kjk =   &K [(j + n*k)*ks1*ks2*ks3];

				for (int i3 = 0; i3 < kk3; ++i3)
				  for (int i2 = 0; i2 < kk2; ++i2)
					for (int i1 = 0; i1 < kk1; ++i1)
					{

						cplxf z = KKtij [i1 + kk1*(i2 + kk2*i3)];

						int j3_0 = max(0, i3-ks3+1), j3_n = min(i3+1, ks3);
						int j2_0 = max(0, i2-ks2+1), j2_n = min(i2+1, ks2);
						int j1_0 = max(0, i1-ks1+1), j1_n = min(i1+1, ks1);

						for (int j3 = j3_0; j3 < j3_n; ++j3)
						  for (int j2 = j2_0; j2 < j2_n; ++j2)
							for (int j1 = j1_0; j1 < j1_n; ++j1)
							{
								z += Kik [(i1-j1) + ks1*
									      ((i2-j2) + ks2*
									       (i3-j3))]
								  * conj (Kjk [(ks1-1-j1) + ks1*
										      ((ks2-1-j2) + ks2*
										       (ks3-1-j3))]);
							}

						KKtij [i1 + kk1*(i2 + kk2*i3)] = z;
					}
			}
		}
	}
}

